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Section 1 

Introduction to registration, integration 
and fusion of remotely sensed data 


Section la 

Overview of issues and challenges 


Image Registration 

in the Context of Space Missions 





Image Registration 

in the Context of Space Missions 




Image Registration 

in the Context of Earth Remote Sensing 


Example of Various Spatial and Spectral Characteristics 
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What is Image Registration . . . 


• Definition 

“Exact pixel-to-pixel matching of two different 
images or matching of one image to a map ” 

• Navigation or Model-Based Systematic Correction 

- Orbital, Attitude, Platform/Sensor Geometric Relationship, Sensor 
Characteristics, Earth Model, etc. 

• Image Registration/Feature-Based Precision 
Correction 

— Navigation within a Few Pixels Accuracy 

- Image Registration Using Selected Features (or Control Points) to 
Refine Geo-Location Accuracy 

• Image Registration as a Post-Processing or as a 
Feedback to Navigation Model 


What is Data Integration and Fusion . . . 


• Definition “Data Integration” 

“Data integration means combining data coming from 

different sources and providing users with a unified view of 
these data” 

— It usually does not refer to the extraction of relevant information, 
but rather to combination and/or concatenation of data (e.g. for a 
GIS or data assimilation or within a model) 

• Definition “Data Fusion” 

“Data fusion is a formal framework in which are expressed 
the means and tools for the alliance of data originating 
from different sources. Data fusion aims at obtaining 
information of greater quality; the exact definition of 
'greater quality' will depend upon the application (Wald, 
1999) 


The Role of Image Registration and Fusion 

in the Processing of Remotely Sensed Data 


• Essential for spatial and radiometric calibration of 
multitemporal measurements for creating long-term 
phenomenon tracking data 


• Used for accurate change detection and land cover/ 
land use assessment 

• Basis for extrapolating data throughout several scales 
for multi-scale phenomena (distinguish between 
natural and human-induced) 


Impact of Misregistration 


• (Towsnhend et al, 1992) and (Dai & Khorram, 1998): small error in registration 
may have a large impact on global change measurements accuracy 

• e.g., 1 pixel misregistration error => 50% error in Vegetation Index (NDVI) 
computation (using 250m MODIS data) 



Human-induced land cover changes observed by Landsat-5 in Bolivia in 1984 and 1998 
(Courtesy: Compton J. Tucker and the Landsat Project, NASA Goddard Space Flight Center) 

• Influence of image registration on products validation 

• Impact of misregistration on legal, economic and sociopolitical (e.g., resource 
management), etc. 




Image Registration and Fusion Applications 


• Multimodal registration and fusion, for integrating 
complementary information from multiple sensors 

• Example: LIDAR and optical data for precision landing 

• Multitemporal registration and integration, for change 
detection and Earth resource surveying 

• Example: Landsat Program, Landsat-4/5, Landsat-7 and Landsat-8 

• Viewpoint registration and fusion, for landmark navigation, 
formation flying (sensor web) and planet exploration 

• Example: Landsat-7, EO-1 and MODIS for disaster management 

• Template registration and integration, for content-based 
searching or map updating 

• Example: Finding flux towers in MODIS or EO-1 data 


Challenges in Image Registration 

for Remote Sensing 


• Remote Sensing vs. Medical or Other Imagery (“Big Data”) 

— Variety in the types of sensor data and the conditions of data acquisition 
— Size of the data 
— Lack of a known image model 
— Lack of well-distributed “fiducial points” 

• Navigation Error (or varying “Initial Conditions”) 

— Historical satellites (e.g., Landsat-5 compared to Landsat-7) 

— Following a maneuver (e.g., star tracking) 

• Needs: 

— Sub-pixel accuracy 
- Robustness to recurring use 

— Speed and High-Level of Autonomy (Near- or Near-real time applications, 
e.g., disaster management) 

— On-the-ground or On-Board Processing 


Challenges in Image Integration and Fusion 
for Remote Sensing 


• Same challenges than registration: 

— Various spectral and spatial resolutions 

— Acquired at different times, under different conditions 

— For example: Time series over a large number of years involving multiple 
programs and/or instruments 

• Misregistration error 

• Difficulty in quality assessment because of the lack of reference 
against which one may compare properties of the fused image; 
validation often done by the end-application not at the time of the 
fusion 

— Objective function often missing 
— Measures such as correlation, entropy, etc. 


Challenges with 

Atmospheric and Cloud Interactions 



Baja Peninsula , California; 4 different times of the day (GOES-8) 


Challenges with 

Multitemporal Effects 

Mississippi and Ohio Rivers before & after Flood of Spring 2002 (Terra/MODIS) 



April 25, 2C02 Mayl8,2DD2 


Challenges with 
Relief Effect 


SAR and Landsat-TM Data of Lope Area, Gabon, Africa 





Challenges with 

Relief Effect (2) 

Use Digital Terrain Models (DEMs) => Taking terrain 
into account in matching 







Theoretical vs. Operational Approaches 
for Image Registration 

• Many promising theoretical approaches with good results on 
specific datasets, but no gold standard algorithm for 
operational use 

• Every instrument has a working operational approach, often 
solved the same way (e.g., with Normalized Cross Correlation) 

• Operational Team Requirements: 

— Know models of sensor/platform/etc. 

— Have access to complete data set 

— Have continuing demands/responsibility 

— Are registering same plots of land again and again - can invest effort in 
data preparation 

— Can’t take big risks on not fully proven methods 

• Do not need one magic method - need toolbox of many 
approaches 


Challenges with 

Institutional Challenges 

• Different communities/literature/requirements 

— Photogrammetry 

- Computer vision/image processing 

- Operational teams 

- Remote sensing/Earth scientists/end users 

• Demanding/varying mission requirements 

— Caution in system design, new methods 

• Difficulty to share data or models between instruments 

• Understanding the requirements 

— Know what data is used for 
— Have to fuse many data sets 
— Have access to ancillary data 

- Know cultural and historical data 


Precision Correction in Operational Systems 

Some Examples - Highlights 


• A VHRR: AUTONAV algorithm computes attitude corrections using Maximum Cross- 
Correlation (MCC) method between sequential images 

• GOE S/METE O SAT: CPs and NOAA Shoreline database (GSHHS) used to match edges 
extracted from meteorological images 

• LAND SAT: CP image chips (lm orthorectified) using Gaussian pyramid, automatic Moravec 
window extraction and NCC or Mutual Information 

• MISR: Database of 120 GCPs (each a collection of nine geolocated image patches of a well- 
defined and easily identifiable ground features, from Landsat, terrain-corrected, data) &ray 
casting simulation software 

• MODIS: Biases and trends in the sensor orientation determined from automated control point 
(CP) matching and removed by updating models of the spacecraft and instrument orientation; 
finer CGPs from Landsat TM and ETM aggregated using PSFs and correlated with NCC 

• SEA WIFS: Reference catalog of islands GCPs and matching using spectral classification and 
clustering of data, “nearest neighbor” and pattern matching techniques 

• SPOT: Reference3 D™ using DEM ortho-rectified simulated reference image in focal plane 
geometry, matching of input image to simulated using NCC and resampling into a cartographic 
reference frame 




VEGETATION: Database of CPs from SPOT for VEGETATION 1 and VEGETATION 1 for 
VEGETATION2; Matching by NCC 


Precision Correction in Operational Systems (2) 


• Operational Environment 

- Platform/sensor models integrated 

- Historical data available for statistics/modeling 

- Robustness and consistency over time is a requirement 

• General Characteristics 

- Use database of Ground Control Points (GCP) or Chips 

- Normalized Cross-Correlation (NCC) is the most common similarity measure 

- Digital Elevation Model (DEM) is rarely integrated in the registration process 

- Cloud masking usually integrated 

- Errors in the [0.15-0.5] range 

• Various approaches. No gold standard approach: 

- Create framework to validate new image registration components and 
algorithms 

- For each algorithm, define “region of convergence” and “region of 
divergence” 

- Provide guidance/recommendations for utilization of algorithms and their 
components 

- Provide fast algorithms for real-time/near-real-time and on-board applications 


Section 1 b 

survey of registration, integration and 

fusion methods 


Brief Image Registration Survey 


Image Registration Frameworks 


• Mathematical Framework 

— Il(x,y) and I2(x,y): images or image/map 

- find the mapping (f,g) which transforms It into 12: 

I2(x,y) = g(Il(fx(x,y),fy(x,y)) + n(x,y) 

» f : spatial mapping 

» g: radiometric mapping 

» n: sensor and other imaging noise 

- Spatial Transformations 

- Translation, Rigid, Affine, Projective, Perspective, Polynomial, ... 

- Radiometric Transformations “g” (Resampling) 

- Nearest Neighbor, Bilinear, Cubic Convolution, . . . 

• Algorithmic Framework (Brown, 1992) 

1. Feature Extraction 

2. Feature Matching (Similarity Metrics & Matching Strategy) 

3. Image Resampling (if needed) 


Image Registration Components 


0 Pre-Processing 

• Cloud Detection, Region of Interest Masking, ... 

1 Feature Extraction (“Control Points”) 

• Gray Levels, Salient Points (e.g.. Edges, Edge-like such as Wavelet 
Coefficients, Corners), Lines, Contours, Regions, Scale Invariant 
Feature Transform (SIFT), etc. 

2 Feature Matching 

• Choice of Spatial Transformation (function f: a-priori knowledge) 

• Choice of Search Strategy : 

Global vs Local, Multi-Resolution, Optimization, ... 

• Choice of Similarity Metrics 

L2-Norm, Normalized Cross-Correlation, Mutual Information, 
Hausdorff Distance, ... 

3 Remapping/Resampling (function g: if necessary) 


General Approaches to Registration 


1 . Manual Registration 

2. Correlation-Based Methods 

3. Fourier-Domain and Other-Transform Based 
Approaches 

4. Mutual Information and Distribution-Based 
Approaches 

5. Feature-Point Methods 

6. Contour- and Region-Based Approaches 


Feature Extraction 


• Gray levels 

• Salient points 

— Edge-like, wavelet coefficients (Simoncelli and Freeman 
‘95) 

- Comers (Kearny et al. ‘87, Harris and Stephens ’88, Shi 
and Tomasi ‘94) 

• Lines 

• Contours, regions (Govindu et al. ‘99) 

• Scale invariant feature transform (SIFT), Lowe ‘04 


Similarity Metrics 


Z, 2 -norm: 


— Minimize the sum of squared errors (SSD) over overlapping 

subimage iV _, 

SSD(x, v) = ^ ^ |/i (m, n)-I 2 (m-x,n- y)J 

m=0 n=0 

Cross-correlation 

— Maximize cross-correlation over image overlap 
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Similarity Metrics (2) 


• Mutual information (MI): 


Maximizes the degree of statistical dependence between the 
images 
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where Mis the sum of all histogram entries, i.e., number of 
pixels (in overlapping subimage) 


Similarity Metrics (3) 


• Partial Hausdorff distance (PHD): 

H k (I,J 2 ) = K* ai min^ dist ( Pl ,p 2 ), 

where 1< K< |/, (Huttenlocher et al. ‘93, Mount et al. ‘99) 



Similarity Metrics (4) 


• Discrete Gaussian mismatch (DGM): 
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is similarity measure ranging between 0 and 1 


Transformation Functions 


• Translation-only, rigid 

• Rotation, scale, and translation (RST) 

• Affine (6 degrees of freedom) 
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• Projective/homography (e.g., for perspective effects in image 
mosaicing); 8 parameters 


Transformation Functions (2) 


• Weighted linear transformation; adaptive transformation, 
continuous and smooth, applied to multiview images with 
local geometric differences, and maps an entire image to 
another 

— Interpolating surface is a weighted sum of planar patches, 
each of which passes through a control point and provides 
a desired gradient, i.e., 


for monotonically decreasing weight R : (x, y) = |(x - v f +( r- r f 




and 


A (*> y) = a i ( x ~ x i ) + b i (y ~ y ) + F i 


Matching Strategies 


• Exhaustive search (exponential in dimensionality of space) 

• Fast Fourier transform (FFT) 

• Numerical optimization (e.g., steepest gradient descent wrt 
SSD, NCC, and MI (Thevenaz, Ruttimann, and Unser (TRU) 
‘98; Spall ’92) 

• Robust transformation estimate (e.g., RANSAC, LMS) if 
(most) correspondences are known (via SIFT-like) 

• “Correspondenceless”, e.g., correlation of descriptor 
distribution/feature consensus (Govindu et al. ‘99) 

• Robust feature matching (RFM), e.g., efficient subdivision and 
pruning of transformation space; Huttenlocher et al. ‘93, Mount 
et al. ’99, Netanyahu et al. ‘04 


Matching Strategies (2) 


• Frequency domain-based approach 

— Efficient computation of correlation as inverse of 
F*(u,v)F 2 (u,v) 

— Practical implementation (extension to NCC, masking 
invalid pixels, optimized computation) 

- Finding (small) rotational and scale differences (by 
matching chips) 

— Subpixel registration for translation-only using phase 
estimate (also in case of image aliasing) 

- Rotation and scale estimate by casting to log-polar 
coordinates 


Matching Strategies (3) 


• Matched filtering 

- Maximize SNR (using theory of linear systems) 

- Apply phase-only and symmetric phase-only matched 
filters for translation-only IR 


— Apply Fourier-Mellin transform for rotation and scale 
changes; transform represents these parameters as 
translational shifts in log-polar coordinates of magnitude of 
Fourier spectrum, i.e., first estimate rotation and scale, 
followed by translation estimate 


Phase product = 


(w,v) F 2 (u,v) _ ~j(ui x +vi y ) 



Matching Strategies (4) 


• Numerical optimization 

- Powel’s, Brent’s (1-D), simplex, etc. 


- Steepest descent/ascent variants 

• Standard P^+i = Pa: “ A 8 k 

• Newton-Raphson P^+i = 

• Levenberg-Marquardt p* +1 = P/ - (l^ + ^diagjH* ]) 1 g k 

- Apply to various similarity metrics, e.g., SSD (Eastman and 
Le Moigne ‘01), Mutual Information, etc. 


» Explicit computation of gradient (and Jacobian/Hessian), e.g., 
Thevenaz and Unser ‘00 

» Stochastic approx. (Spall ‘92; Cole-Rhodes et al. ’03) 


Matching Strategies (5) 


• Alignment via local geometric distributions 



Slope angle distributions and their correlation 


Matching Strategies (6) 


• Robust feature matching (RFM) 

- Space of affine transformations: 6-D space 

- Subdivide: Quadtree or kd-tree. Each cell T represents a set 
of transformations; T is active if it may contain/ ; o/w, it is 

killed 

— Uncertainty regions (UR’s): Rectangular approximation to 
the possible images r(a) for all rGT,a 

— Bounds: Compute upper bound (on optimum similarity) by 
sampling a transformation and lower bound by computing 
nearest neighbors to each UR 

- Prune: If lower bound exceeds best upper bound, then kill the 
cell; o/w, split it 


Matching Strategies (7) 


• Computational efficiency 

- “Culling” feature points via, e.g., condition theory 

— Efficient numerical or discrete algorithmic 
procedures 

- Hierarchical pyramid-like (wavelet) decomposition 

- Use landmark chip database (instead of a large 
scene) or alternatively, extract automatically 
corresponding regions of interest using 
mathematical morphology (Plaza et al ‘07) 


Brief Image Fusion Survey 


Methods of Image Fusion 


• Prerequisite to Image Fusion: Very accurate registration 

• Various Approaches [Blum and Liu, 2005]: 

— Multiscale-Decomposition-Based Fusion Methods 
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1. Multiscale representations: Pyramid Transform (e.g., Laplacian), Discrete Wavelet 
Transform (DWT) and Discrete Wavelet Frame (DWF) 

2. Activity-level measurement 

3. Coefficient Grouping Method 

4. Coefficient Combining Method 

5. Consistency Verification 






Methods of Image Fusion (2) 


- NonMultiscale-Decomposition-Based Methods 

• Pixel-Level Weighted Averaging 

- e.g., PCA 

• Non-Linear Method 

- e.g., separate images into low-pass (LP) and high-pass (HP); then 
modify and fuse LP and HP 

• Estimation Theory Based Methods 

- e.g. using Maximum A Priori (MAP) and Maximum Likelihood 
(ML) estimates or Markov Random Field (MRF) distributions 

• Color Composite Fusion 

- Combines input images in color space => false color representation 

- Often used with another method, e.g., PCA or Neural Network 

• Artificial Neural Networks 

- e.g., trained to superimpose objects of interest on background 


Performance Evaluation of Fusion Algorithms 


Objective Evaluation using a Reference Image 

• Root Mean Square Error ( RMSE) 


RMSE = 


i 


i NM 


NM Z 
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» Where R is the reference image, F the fused image, and NxM 
the size of the image 

Correlation 

Peak Signal to Noise Ratio ( PSNR ) 


PSNR = 10 log 
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» Where L is the number of gray levels in the image 


Performance Evaluation of Fusion Algorithms (2) 


• Objective Evaluation using a Reference Image 
(cont.) 

• Mutual Information (MI) 

• Universal Quality Index (QI) (Wang and Bovik, 2002) 

4a xy 

Q (a 2 x+ a 2 y )[(x) 2 + (y) 2 ] 

» Where x and y represent the average of all pixels in the 

reference and the fused images, respectively, 

2 2 

» o x and o y represent the respective variances and o ^ represents 
the co-variance between the 2 images 

- Other indexes derived from QI, (e.g., (Piella et al, 2003) where 
index based on the local saliencies of the input images) 


Performance Evaluation of Fusion Algorithms (3) 


• Objective Evaluation without a Reference Image 

• Standard Deviation (SD) 

• Entropy (H) 

• Overall Cross Entropy of the source images X, Y and the 
fused image F 


CE(XJ-F) 
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CE(X-F) = ^h x (i)\og 2 


i = 0 


%(!) ' 
\ (0 ) 


» Where h is the normalized histogram of the image 


References 


• • • 


• L. Wald, “Some terms of reference in data fusion,” IEEE Transactions on Geoscience and Remote Sensing, 
Vol. 37, No. 3, 1999, pp. 1190-1193. Also at http://www.data-fusion.org/terms_of_reference 

• J.R. Townshend, C.O. Justice, C. Gurney, and J. McManus, “The impact of misregistration on change 
detection,” IEEE Transactions on Geoscience and Remote Sensing, Vol. 30, No. 5, 1992, pp. 1054-1060. 

• X. Dai and S. Khorram, “The effects of image misregistration on the accuracy of remotely sensed change 
detection,” IEEE Transactions on Geoscience and Remote Sensing, Vol. 36, No. 5, 1998, pp. 1566-1577. 

• L.G. Brown, “A survey of image registration techniques,” ACM Computing Surveys, Vol. 24, No. 4, 1992, 
pp. 325-376. 

• J. Le Moigne, N.S. Netanyahu, and R.D. Eastman, Image Registration for Remote Sensing, Cambridge 
University Press, 2011. 

• E.P. Simoncelli and W.T. Freeman. “77ze Steerable Pyramid: A Flexible Architecture for Multi-Scale 
Derivative Computation, ” IEEE Second Int'l Conf on Image Processing, October 1995. 

• J.K. Kearney, W.B. Thompson and D.L. Boley, “Optical flow estimation: An error analysis of gradient-based 
methods with local optimazation,” IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. 9 , 
1987, pp. 229-244. 

• C. Harris and M. Stephens, “A combined comer and edge detector,” Proceedings of the 4 th Alvey Vision 
Conference, 1988, pp. 147-151. 

• J. Shi and C. Tomasi, “Good features to track,” 9 th IEEE Conference on Computer Vision and Pattern 
Recognition, Springer. 

• V. Govindu and C. Shekhar. Alignment Using Distributions of Local Geometric Properties, IEEE 
Transactions on Pattern Analysis and Machine Intelligence, PAMI, October 1999. 

• D.G. Lowe, “Distinctive Image Featyres from Scale-Invariant Keypoints,” International Journal of Computer 
Vision, Vol. 60, No. 2, 2004, pp. 91-110. 


References 


• • • 


• D.P. Huttenlocher, G.A. Klanderman and WJ. Rucklidge, “Comparing images using the Hausdorff distance,” 
IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. 15 , 1993, pp. 850-863. 

• D.M. Mount, N.S. Netanyahu, and J. Le Moigne, “Efficient algorithms for robust point pattern matching’, 
Pattern Recognition, Vol. 32 , 1999, pp. 17-38. 

• N.S. Netanyahu, J. Le Moigne and J.G. Masek, “Georegistration of Landsat data via robust matching of 
multiresolution features,” IEEE Transactions on Geoscience and Remote Sensing, Vol. 42 , 2004, pp. 
1586-1600. 

• P. Thevenaz, U. Ruttimann and M. Unser, “A pyramid approach to subpixel registration based on intensity,” 
IEEE Transactions on Image Processing, Vol. 7, No. 1, 1998, 27-41. 

• J.C. Spall, “Multivariate stochastic approximation using a simultaneous perturbation gradient 
approximation,” IEEE Transactions on Automatic Control, Vol. 37, No. 3, 1992, pp. 332-341. 

• A. Cole-Rhodes, K. Johnson, J. Le Moigne, and I. Zavorin, “Multiresolution registration of remote sensing 
imagery by optimization of mutual information using a stochastic gradient,” IEEE Transactions on Image 
Processing, Vol. 12, No. 12, 2003, pp. 1495-1511. 

• A. Plaza, J. Le Moigne, and N.S. Netanyahu, “Parallel Morphological Feature Extraction for Automatic 
Registration of Remotely Sensed Images,” 2007 IEEE International Geoscience and Remote Sensing 
Symposium, IGARSS'07, Madrid, Spain, July 2007. 

• R.S. Blum and Z. Liu (eds), “Multi-Sensor Image Fusion and its Applications,” Taylor & Francis, 2005. 

• Z. Wang and A.C. Bovik, “A Universal Image Quality Index,” IEEE Signal Processing Letters, Vol. 9, No. 3, 
2002, pp. 81-84. 

• G. Piella and H. Heijmans, “A New Quality Metric for Image Fusion,” Proc. International Conference on 
Image Processing (ICIP), Sept 2003, Vol.2, pp. Ill- 173- 176. 



Section 2c 

Wavelets and Redundant Representations 

for Image Registration 


Why Wavelets and Redundant Representations for 
Image Registration 

• Representations such as the multiresolution wavelets utilized for the new 
compression standard JPEG-2000 can bring multisensor/multiresolution 
data to the same spatial resolution without losing significant information 
and without blurring the higher resolution data. 

• At the lower resolutions, the process preserves important global features 
such as rivers, lakes and mountain ridges as well as roads and other man- 
made structures, while at the same time eliminating weak higher 
resolution features often considered as "noise or "spurious pixels." Of 
course, some of the important finer features will be lost at the lowest 
resolutions but, if needed, they can be retrieved in the iterative higher 
levels of decomposition. 

• The multiresolution iterative search focuses progressively toward the 
final transformation with a decreasing search interval and an increasing 
accuracy at each iteration. Hence, this type of strategy achieves higher 
accuracies with higher speeds than a full search at the full resolution. 

• Wavelet decomposition and multiresolution iterative search are very 
well-suited for fme-grained parallelization, thus speeding up the 
computations even more. 


Wavelets for Image Registration 


• 1994: First results on the utilization of orthogonal 
Daubechies wavelets for image registration 





Choice of 


Best 


TraDsforiuatioD 



Figure 1 
Original Image 



Figure 2 Figure 3 

Wavelet Coefficients Corresponding Wavelet Coefficients Correspond 

to Figure 1 to Figure 1 rotated 44 degrees 



Rotation and Translation Invariance Issues 


• Orthogonal wavelets, feature information changes within 
or across subbands with subsampling => Study for Shift 
Sensitivity: 

- low-pass subband relatively insensitive to translation, if 
features are twice the size of wavelet filters 

— high-pass subband more sensitive but can still be used. 
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Correlate Wavelets of Two Pulses 


Translation Sensitivity - Low- Pass Level 3 


Correlation Minima vs Pulse width 
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Translation Sensitivity - High-Pass Level 3 


Correlation Minima vs Pulse width 
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Correlation Averages vs Pulse width 



Rotation- and Translation-Invariant 
Representations 


• Spline Wavelets [Battle & Lemarie; Unser et al] 


V" = fg" (x) = ^ c t ( k)cp n (2"' x - k) , x G 9t, c ; G l 2 } with scaling function (p n (x) = ^ p(k)P\x-k) 

k =- oo 

p arbitrary invertible convolution operator or filter, 
and P n (x) is a 5-spline of order n (can be constructed 
by repeated convolution of 5-Spline of order 0) 


Example of B-Spline Scaling Function and Associated Wavelet 


k =- oo 


n=3 (cubic) 



Original 
Image or 
Output of LI 


•© 


Simoncelli et al 

- Relax critical sampling condition of wavelet transforms 

- Provides an overcomplete representation by 4k/3 





Framework for Evaluation of 

Image Registration Components 











Framework for Evaluation of 

Image Registration Components 


Features 


Similarity 

Measure 


Strategy 
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Algorithm Testing Using 
Synthetic Data 


• • • 



Transformation of Starting Scene: 

• Scales in [0.8, 1.2] (step = 0.05) 

• Translations in [0,20] pixels (step = 0.5) 

• Rotations in [0,20] degrees (step = 0.5) 

• Gaussian noise in [-15, 20] dB (step = 1) 

• Radiometric Transformation (PSF 
constructed from black 512x512 image 
with 5x5 white center) 



Synthetic Image Examples (Original; Warp & Noise; Warp & PSF) 





Results Marquart Levenberg Optimization + L2 Similarity 
Spline Wavelets and Simoncelli Features 
Varying Noise - With or Without PSF [Zavorin et Le Moigne, 2005] 



Warping + Noise 


TFIU-3fe«C THU^-SimB TflU-atmL 



Warping + PSF 


The shaded green regions in each plot correspond to those 
points in the parameter space for which the resulting global 
RMS registration error is larger than 1.0 


TRU-S04C TRo-SimB TRu-£hnL 



Warping + PSF + Noise 








Sensitivity to Initial Guess [Le Moigne et al, 2004] 


THU — SpIC E(p ) p- 1 


THU-sime i 


TRU— SifflL E(p t > P- 1 



Sensitivity of TRU Algorithms to Initial Guess 
(with PSF) 


Mt-SplG >■ 1 


MI-SimB >■ 1 


Ml— Sim l_ 


Overall Main Findings 


Simoncelli-based methods outperform 
Spline pyramid-based methods 
Optimization based on Mutual Information 
does not perform better than L2-Norm 
Simoncelli/Low-Pass better than Simoncelli/ 
Band-Pass for Low Noise and Same 
Radiometry and for Initial Guess Sensitivity 
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Sensitivity of TRU/MI to initial Guess 








Registration using a Chip Database 





1. Find Chips that correspond to 
the Incoming Scene 

2. For Each Chip, Extract Window 
from input scene using UTM 
coordinates 

3. Eliminate Windows with 

insufficient information 

4. Smooth and Normalize gray 
values of both Chip and 
Window using a Median Filter 

5. Register each ( Chip, Window) 
Pair using a wavelet-based 
automatic registration: get a 
local transformation for each 
pair 

6. Eliminate Outliers 

7. Compute Global Rigid 

Transformation from all local 
ones 

8. Compute Correct UTM of 4 
Scene Corners of input scene 

9. If desired, Resample the input 
scene according to the global 
transformation 


Algorithm Testing Using . . . 

Landsat-TM Multitemporal Data [Netanyahu et al, 2004] 
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Algorithm Testing Using . . . 

Landsat-TM Multitemporal Data (2) 


• Landsat-5 and -7 Multi-Temporal Data [Netanyahu et al, 2004]: 
Chips and Corresponding Windows 




One chip and 4 Corresponding 
Windows Extracted from 4 
Multi-Temporal Landsat Imagery 


Algorithm Testing Using . . . 

Landsat-TM Multitemporal Data (3) 


Chip # 

84240 

87136 

96193 

97275 

#1 - Rot: 

0 

0 

1 

0 

TX 

8 

18 

12 

21 

TY 

-39 

-25 

-92 

-28 

Distance 

0.00 

1.00 

1.80 

0.00 

#2 - Rot: 

0 

-0.5 

1 

0 

TX 

8 

11 

-66 

22 

TY 

-41 

-41 

4 

-30 

Distance 

0.00 

0.62 

1.93 

0.00 

#3 - Rot: 

0 

-0.5 

-0.3 

0 

TX 

8 

11 

10.84 

21 

TY 

-40 

-41 

-96 

-29 

Distance 

1.00 

0.89 

0.52 

0.00 

#4 - Rot: 

0 

0.8 

0 

0 

TX 

8 

12.34 

11 

21 

TY 

-41 

-38.12 

-94 

-28 

Distance 

1.00 

0.96 

0.00 

0.00 

#5 - Rot: 

0.4 

0.7 

0 

0 

TX 

7.8 

10.4 

-38 

20 

TY 

-40.86 

-40.34 

52 

-31 

Distance 

0.94 

0.93 

2.24 

0.00 

#6 - Rot: 

0 

-0.8 

0.3 

0 

TX 

6 

10.61 

11.5 

22 

TY 

-41 

-41.94 

-99 

-33 

Distance 

1.00 

0.78 

0.70 

0.00 

#7 - Rot: 

0.5 

0 

0 

0 

TX 

9 

12 

12 

22 

TY 

-40 

-38 

-94 

-29 

Distance 

0.56 

0.00 

0.00 

0.00 


Transf. 

84240 

87136 

96193 

97275 

Rotation 

0.013 

0.003 

-0.042 

-0.143 

Transl-x 

7.18 

11.43 

12.61 

21.20 

Transl-y 

-41.12 

-40.49 

-95.38 

-28.85 



Global Registration for 4 Scenes 


Transf. 

84240 

87136 

96193 

97275 

Rotation 

0.00 

0.00 

0.00 

0.00 

Transl-x 

7.18 

10.55 

9.48 

20.97 

Transl-y 

-40.06 

-39.16 

-95.16 

-28.97 


Manual Registration for 4 Scenes 


• This Chesapeake Bay Example: 

Global Accuracy Error ~ 0.82 pixel 

• Other Virginia Scenes: 

Global Accuracy Error ~ 0.31 pixel 


Algorithm Testing Using . . . 

EO-1 and Global Land Survey (GLS) Maps 



Algorithm Testing Using . . . 





Algorithm Testing Using . . . 
Multisensor Data 


ETM+ 

• Multi-Sensor Data 

- EOS Validation Core Sites 

- IKONOS/Landsat-7/MODIS/ 

SeaWiFS 

• Red and NIR bands for each sensor 

• Spatial resolutions: IKONOS: 4m; 

ETM+: 30m; MODIS: 500m; 

SeaWiFS: 1000m ETM/IKONOS - Coastal VA Data 


— 4 different sites: 

• Coastal Area: VA, Coast Reserve 
Area, October 200 1 

• Agriculture Area: Konza Prairie in 
State of Kansas, July to August 2001 

• Mountainous Area: Cascades Site, 

September 2000 

• Urban Area: USDA Site, Greenbelt, 

MD, May 200 1 ETM/IKONOS - Agricultural 

Konza Data 



ETM+ 


IKONOS 



Algorithm Testing Using . . . 

EOS Validation Core Sites (2) 


Pair to Register 

Method 1 (GC) 
Rotationj Translation 

Method 2 (GGD) 
Rotationj Translation 

Methoc 

Rotation 

l 3 (WCE) 
Translation 

Method 

Rotation 

4 (WMIE) 
Translation 

Methoc 

Rotation 

l 5 (WHR) 
Translation 

(1) etm_nir_3 1.25. power / 
etm_red_31. 25. extract 


Rot 

ation = 0 , ' 

Translation = (0,0 

>) computed 

by all methods, 

using sever 

i sub-windows p 

airs 


(2) iko_nir_3. 91. power / 
etm_nir_3 1.25. extract 

- 

(2,1) 

0.0001 

(1.9871,-0.0564) 

0 

(2,0) 

0 

(2,0) 

0 

(0,0) 

(3) iko_red_3.91.power / 
etm_red_31. 25. extract 

- 

(2,1) 

-0.0015 

(1.7233,0.2761) 

0 

(2,0) 

0 

(2,0) 

0 

(0,0) 

(4) etm_nir_31.25.power / 
modis_day249_cc_nir .extract 

- 

(-2,-4) 

0.0033 

(-1,7752,-3.9238) 

0 

(-2,-4) 

0 

(-2,-4) 

0 

(-3,-3 .5) 

(5) etm_red_31.25.power / 
modis_day249_cc_red. extract 

- 

(-2,-4) 

0.0016 

(-1.9665,-3.9038) 

0 

(-2,-4) 

0 

(-2,-4) 

0 

(-2,-3 .5) 

(6) modis_day249_cc_nir .power / 
sea wifs_day256_to249_nir. extract 

- 

(-9,0) 

0.0032 

(-8.1700,0.2651) 

0 

(-8,0) 

0 

(-9,0) 

0.5 

(-62) 

(7) modis_day249_cc_red.power / 
seawifs_day256_to249_red. extract 

- 

(-9,0) 

0.0104 

(-7.6099,0.5721) 

0 

(-8,0) 

0 

(-8,0) 

0.25 

(-7,1) 


• GC: Gray Levels + Fast Fourier Correlation 

• GGDL Gray Levels + Gradient Descent 

• WCE: Wavelets + Correlation 

• WMIE: Wavelets + Mutual Information 

• WHR: Wavelets + Hausdorff + Robust Feature Matching 


Algorithm Testing Using . . . 

EOS Validation Core Sites (3) 


Image Name 

Computed X 

Computed Y 

Comes from 
Registered Pair 

IKONOS red 

0 

0 

(Starting Point) 

IKONOS nir 

-0.2500 

-0.2500 

IKO red to ETM red 
and ETM red to IKO nir 

1 IKONOS nir 

-0.2500 

-0.3125 

IKO red to ETN nir 
and ETM nir to IKO nir 


Table 3 - Self-Consistency Study of the Normalized Correlation Results 


Image Name 

Computed X 

Computed Y 

Comes from 
Registered Pair 

IKONOS red 

0 

0 

(Starting Point) 

IKONOS nir 

0.2500 

0.0000 

IKO red to ETM red 
and ETM red to IKO nir 

IKONOS nir 

0.1250 

-0.1250 

IKO red to ETN nir 
and ETM nir to IKO nir 


Table 4 - Self-Consistency Study of the Mutual Information Results 


These consistency studies show between 0.125 and 0.25 pixel errors using circular 

registrations of IKONOS NIR and Red data 
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Section 2f 

Fusion using cokriging 


Kriging and Cokriging 


• • • 


• Kriging and cokriging are geostatistical techniques 
used for interpolation purposes, originally used in 
geo-statistics, mining, and petroleum engineering 
applications 

• Both methods are generalized forms of univariate and 
multivariate linear regression models: 

— They are linear-weighted averaging methods, similar to 
other interpolation methods 

- Other methods use weights based on the distance of each 
control point (sample value) from the target location; 
controls points closer to the target receive the larger 
weights => not necessarily true if the data exhibit strong 
anisotropy 


Kriging and Cokriging (2) 


• Kriging and cokriging 

- Pioneered by Danie Krige, 1951; formalized by Georges 
Matheron in the 1 960s 

— Ability to capture anisotropy of the underlying variables 
through the spatial covariance model: 

• Distant control points along the axis of maximum correlation should 
have greater influence on the interpolated value 

• Weights depend not only on distance, but also on the direction and 
orientation of the neighboring data to the unsampled 

• Generalized version of kriging (B.L. U.E): 

— Best: aims to minimize variance of the errors 

— Linear: estimates are weighted linear combination of the 
available data 

- Unbiased: tries ta have mean residual, or error, equal to zero. 
— Estimator: k = V Wj.v 


Cokriging 


• • • 


Interpolation using more that one type 
of variable to estimate an unknown 
value at a particular location. 

Estimation error: 


Goal of cokriging is to minimize 
variance of error subject to some 
constraints (to ensure unbiasedness of 
our estimate, here “ordinary cokriging”) 


n m 

ii 0 = '^a i .ii l + Yb J .v j 

R = U 0 -U 0 =w'Z, 
w‘ =(a l ,...,a n ,b l ,...,b m ,~ 1 ) 
Z‘ = (U l ,...,U n ,V l ,...,V m ,U 0 ) 


Var(R) = W C z w 


n m 



Cokriging (2) 


• Variance to minimize with 2 constraints: 


Var(R) = w‘C z w 

VI VI Wl Wl 

- 2 2 i.YV b : b ,Cov<,yy > 

i j i j 


VI TVl VI 

+2 22 a Jj^oviUy^ - 2 ^ a,Cov(U,U 0 ) 

i j i 

m 

- 2 Y bjCov(VjU 0 ) + Cov(U 0 U 0 ). 

j 


Var(R) = w‘C z w + 2fi x 





Where ml and m2 are 2 Lagrange multipliers 


Cokriging (3) 


• The next step is taking partial derivatives of the above equation 
with respect to all n + m cokriging variables and the two 
Lagrange multipliers and setting them to zero. Then, we have 
the following n+m+2 equations to solve: 


n m 


J a l Cov(U i U J ) + J b l Cov(V i U j )+«, = Cov(U l) U J ) , for (j = 1 . . . n) 


1=1 


i = 1 


2 a l Cov(Uy ] ) + 2 b l Cov(Vy j ) + ft 2 = Cov(U 0 Vj), for (j = 1 ...m) 


i = 1 


/=1 



i=l 



Experiments Cokriging - Spectral Dimension 


• Spectral dimension 

- ALI: 9 multispectral bands 

- Hyperion: 220 hyperspectral bands 

- Increase spectral resolution of ALI where needed using Hyperion 
information 

• One variable only 

• Software used 

— U CL-FAO Agromet proj ect 

(http : / / www. aigeostats . or g/ software/Geo stats_software/ 
agromet.htm) 


C++ 


reflectance 


Results Cokriging - Spectral Dimension 


Original 1 pixel plot for ALI and Hyperion 



reflectance 


Results Cokriging - Spectral Dimension (2) 



Fusion results on one pixel using cokriging by creating one band/value 
in center of each wavelength interval where ALI data is missing. 


wavelength (nm) 


reflectance 


Results Cokriging - Spectral Dimension (3) 


Fusion results on one pixel using cokriging by estimating up to 3 values 
in each wavelength interval where ALI data is missing. 



i I 

— Hyperion 
-0- ALI 
-X- Fused ALI 
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400 1600 1800 2000 2200 
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reflectance 


Results Cokriging - Spectral Dimension (4) 


Fusion results on one pixefusing cokriging by estimating values at all Hyperion centers 

in each wavelength interval where ALI data is missing. 



Experiments Cokriging - 

Fusion of Panchromatic and Multispectral 

• Landsat 7 ETM data sets 

— Datasets provided by IEEE Data Fusion Committee 

- 8 Bands: 

• 7 multispectral at 30 m spatial resolution 

• 1 panchromatic at 1 5 m spatial resolution 

• Objective 

- Pan sharpening of selected multispectral bands 

• Comparative Methods 

- Cokriging 

- PCA 

- Wavelets 


Experiments Cokriging - 

Fusion of Panchromatic and Multispectral (2) 


Landsat 7 Multispectral 
Bands 2, 3, and 4 



Band 

Spatial 

(meters) 

Resolution 

Spectral 

(m) 

1 

30 

0.45-0.52 

2 

30 

0.53-0.61 

3 

30 

0.63-0.69 

4 

30 

0.78-0.90 

5 

30 

1.55-1.75 

6 

30 

10.4-12.5 

7 

30 

2.09-2.35 

8 (PAN) 

15 

0.52-0.90 


Landsat Panchromatic Band 8 



Methods: Cokriging - 

Fusion of Panchromatic and Multispectral 


Spectral Resolution 


Pan + MS-2 
Pan + MS-3 
Pan + MS-4 



MS-1 


MS-2 


MS-3 


MS-4 


MS-5 


MS-7 


i=> fused_b2 
i — > fused_b3 
i=> fused b4 


Spectral Resolution 
1 pixel of MS band 


MS-Value 


Pan Value 

Pan Value 

1 

2 

Pan Value 

Pan Value 

3 

4 

xl yl 

pi ? 

x2 y2 

p2 ? 

x3 y3 

p3 ms 

x4 y4 

p4 ? 


(using nearest 
Neighborhs, e.g. 5x5) 


Methods: Principal Components Analysis - 

Fusion of Panchromatic and Multispectral 


Input Raw Data Principal Components 


Fusion 


Fused Bands 


MS-1 


MS-2 


MS-3 


MS-4 


MS-5 


MS-7 


PC Transform 
■ > 


PAN 


Its mean and standard deviation. 


t 



S- 

PAN 


PC-1 

Replace PC-1 
with stretched 
PAN band 

S- 

PAN 

Inverse 

PC Transform 

FI 

PC-2 

F2 

PC-2 

PC-3 

F3 

PC-3 


> 

PC-4 

> 


PC-4 

1 v* 

1 1 

F4 

PC-5 


PC-5 


F5 

PC-7 

* m mam ^m m u 

PC -7 


F7 


Methods: Wavelet-Based Fusion - 

Fusion of Panchromatic and Multispectral 


High Spatial 
Resolution 
Data: PAN 


• • • 







Low Spatial 
Resolution 
Data: MS 



WWW 

♦ ♦♦ 


Decomposition 




FUSED 

DATA 





Improved 

Resolutions 


Reconstruction 





Evaluation Fusion Results - 

Fusion of Panchromatic and Multispectral 

• Past Quality Metrics 

- Piella, etc. 

— Gray level only 

- No support for multi-spectral image 

• Objective 

- Improved Classification 

• Performed k-means with k=7, max iterations 1 5 (for PCA and wavelets) 

- Needs ground truth 

• Similarities 

— Spatial Quality: Entropy 
— Spectral quality: correlation 

• Differences 

— Added Texture 

— Co-occurrence matrix for statistical texture properties (Haralick,1973) 
— Variance image 


Results Cokriging - 

Fusion of Panchromatic and Multispectral 



Landsat Pan-sharpened MS bands 2, 3, and 4 
Through Cokriging with Pan band 8 



Results PCA and Wavelets - 

Fusion of Panchromatic and Multispectral 



Landsat Pan-sharpened MS bands 2, 3, and 4 

with Pan band 8 



Results Evaluation though Correlation - 

Fusion of Panchromatic and Multispectral 


Wavelet 

PCA 

Cokriging 

• Fuse pair of bands at a time 

• Require various spatial 
resolution to differ by power 
of 2 

• Fuse multiple bands at a 
time 

• Require the same number 
of pixels (same image 
dimension) for all bands 

• Fuse multiple bands at a time 

• Fuse data with different 
spatial and spectral 
resolutions 

• Fuse data of different natures 

• Can handle scattered data 


CORRELATION OF FUSED BANDS WITH MS INPUT BANDS 


Bands 

Wavelet 

PCA 

Cokriging 

f2, b2 

0.82 

0.99 

0.91 

f3, b3 

0.84 

0.99 

0.93 

f4, b4 

0.92 

0.75 

0.93 

Average 

0.86 

0.91 

0.92 


Results Evaluation though Entropy - 

Fusion of Panchromatic and Multispectral 


ENTROPY OF MS AND FUSED BANDS 


Original 

Bands 


Fused 

Bands 

Wavelet 

PCA 

Cokriging 

b2 

2.68 

f2 

3.12 

2.69 

3.23 

b3 

3.01 

o 

3.28 

3.72 

3.64 

1)4 

3.44 

f4 

3.93 

5.21 

4.90 

Average 

3.04 


3.44 

3.87 

3.92 


Results Evaluation though Co-Occurrence Matrix - 
Fusion of Panchromatic and Multispectral 


MEAN ENTROPY OF ENTROPY IMAGES OBTAINED THROUGH 

CO-OCCURRENCE MATRICES 


Original 

Bands 


Fused 

Bands 

Wavelet 

PCA 

Cokriging 

b2 

1.37 

f2 

1.37 

1.37 

1.44 

b3 

1.42 

f3 

1.45 

1.49 

1.45 

b4 

1.77 

f4 

1.78 

2.02 

1.96 

Average 

1.52 


1.53 

1.63 

1.62 
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